#include <string.h>
#include <math.h>
#include "phg.h"
static const FLOAT *
RT1_bas(DOF *dof, SIMPLEX *e, int no0, int no1, const FLOAT *lambda)
/* evaluation of basis functions */
{
    static FLOAT values[NFace][Dim];
    FLOAT (*p)[Dim] = values;
    FLOAT L0 = lambda[0], L1 = lambda[1], L2 = lambda[2], L3 = lambda[3];
    INT i, a, b, c;
    const FLOAT *J, *g1, *g2, *g3, *g4;

    J = phgGeomGetJacobian(dof->g, e);
#if 0
    FLOAT *f1, *f2, *f3, *f4;
    f1 = phgAlloc(3 * sizeof(FLOAT));
    f2 = phgAlloc(3 * sizeof(FLOAT));
    f3 = phgAlloc(3 * sizeof(FLOAT));
    f4 = phgAlloc(3 * sizeof(FLOAT));
#endif
    FLOAT f1[3], f2[3], f3[3], f4[3];

    if (no1 <= 0)
	no1 = NFace;
    assert(no0 < no1);

    switch (no0) {
	case 0:
	    a = e->verts[1];
	    b = e->verts[2];
	    c = e->verts[3];
	    if ((a < b) && (b < c)) {
		g2 = J + 4;
		g3 = J + 8;
		g4 = J + 12;
		*f2 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f2 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g2 + 2) - *(g4 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g2) - *(g4) * *(g2 + 2);
		*(f3 + 2) = *(g4) * *(g2 + 1) - *(g4 + 1) * *(g2);
		*f4 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f4 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L1 * *(f2 + i) + L2 * *(f3 + i) + L3 * *(f4 + i);
		}
	    }
	    if ((a < c) && (c < b)) {
		g2 = J + 4;
		g3 = J + 12;
		g4 = J + 8;
		*f2 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f2 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g2 + 2) - *(g4 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g2) - *(g4) * *(g2 + 2);
		*(f3 + 2) = *(g4) * *(g2 + 1) - *(g4 + 1) * *(g2);
		*f4 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f4 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L1 * *(f2 + i) + L3 * *(f3 + i) + L2 * *(f4 + i);
		}
	    }
	    if ((b < c) && (c < a)) {
		g2 = J + 8;
		g3 = J + 12;
		g4 = J + 4;
		*f2 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f2 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g2 + 2) - *(g4 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g2) - *(g4) * *(g2 + 2);
		*(f3 + 2) = *(g4) * *(g2 + 1) - *(g4 + 1) * *(g2);
		*f4 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f4 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L2 * *(f2 + i) + L3 * *(f3 + i) + L1 * *(f4 + i);
		}
	    }
	    if ((b < a) && (a < c)) {
		g2 = J + 8;
		g3 = J + 4;
		g4 = J + 12;
		*f2 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f2 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g2 + 2) - *(g4 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g2) - *(g4) * *(g2 + 2);
		*(f3 + 2) = *(g4) * *(g2 + 1) - *(g4 + 1) * *(g2);
		*f4 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f4 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L2 * *(f2 + i) + L1 * *(f3 + i) + L3 * *(f4 + i);
		}
	    }
	    if ((c < a) && (a < b)) {
		g2 = J + 12;
		g3 = J + 4;
		g4 = J + 8;
		*f2 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f2 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g2 + 2) - *(g4 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g2) - *(g4) * *(g2 + 2);
		*(f3 + 2) = *(g4) * *(g2 + 1) - *(g4 + 1) * *(g2);
		*f4 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f4 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L3 * *(f2 + i) + L1 * *(f3 + i) + L2 * *(f4 + i);
		}
	    }
	    if ((c < b) && (b < a)) {
		g2 = J + 12;
		g3 = J + 8;
		g4 = J + 4;
		*f2 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f2 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g2 + 2) - *(g4 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g2) - *(g4) * *(g2 + 2);
		*(f3 + 2) = *(g4) * *(g2 + 1) - *(g4 + 1) * *(g2);
		*f4 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f4 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L3 * *(f2 + i) + L2 * *(f3 + i) + L1 * *(f4 + i);
		}
	    }
	    for (i = 0; i < Dim; i++)
		(*p)[i] = (*p)[i] * 2;
	    if (no1 == 1)
		return (FLOAT *)values;
	    p++;

	case 1:
	    a = e->verts[0];
	    b = e->verts[2];
	    c = e->verts[3];
	    if ((a < b) && (b < c)) {
		g1 = J;
		g3 = J + 8;
		g4 = J + 12;
		*f1 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f1 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f3 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g3 + 2) - *(g1 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g3) - *(g1) * *(g3 + 2);
		*(f4 + 2) = *(g1) * *(g3 + 1) - *(g1 + 1) * *(g3);
		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L0 * *(f1 + i) + L2 * *(f3 + i) + L3 * *(f4 + i);
		}
	    }
	    if ((a < c) && (c < b)) {
		g1 = J;
		g3 = J + 12;
		g4 = J + 8;
		*f1 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f1 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f3 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g3 + 2) - *(g1 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g3) - *(g1) * *(g3 + 2);
		*(f4 + 2) = *(g1) * *(g3 + 1) - *(g1 + 1) * *(g3);
		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L0 * *(f1 + i) + L3 * *(f3 + i) + L2 * *(f4 + i);
		}
	    }
	    if ((b < c) && (c < a)) {
		g1 = J + 8;
		g3 = J + 12;
		g4 = J;
		*f1 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f1 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f3 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g3 + 2) - *(g1 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g3) - *(g1) * *(g3 + 2);
		*(f4 + 2) = *(g1) * *(g3 + 1) - *(g1 + 1) * *(g3);
		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L2 * *(f1 + i) + L3 * *(f3 + i) + L0 * *(f4 + i);
		}
	    }
	    if ((b < a) && (a < c)) {
		g1 = J + 8;
		g3 = J;
		g4 = J + 12;
		*f1 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f1 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f3 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g3 + 2) - *(g1 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g3) - *(g1) * *(g3 + 2);
		*(f4 + 2) = *(g1) * *(g3 + 1) - *(g1 + 1) * *(g3);
		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L2 * *(f1 + i) + L0 * *(f3 + i) + L3 * *(f4 + i);
		}
	    }
	    if ((c < b) && (b < a)) {
		g1 = J + 12;
		g3 = J + 8;
		g4 = J;
		*f1 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f1 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f3 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g3 + 2) - *(g1 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g3) - *(g1) * *(g3 + 2);
		*(f4 + 2) = *(g1) * *(g3 + 1) - *(g1 + 1) * *(g3);
		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L3 * *(f1 + i) + L2 * *(f3 + i) + L0 * *(f4 + i);
		}
	    }
	    if ((c < a) && (a < b)) {
		g1 = J + 12;
		g3 = J;
		g4 = J + 8;
		*f1 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f1 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f3 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g3 + 2) - *(g1 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g3) - *(g1) * *(g3 + 2);
		*(f4 + 2) = *(g1) * *(g3 + 1) - *(g1 + 1) * *(g3);
		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L3 * *(f1 + i) + L0 * *(f3 + i) + L2 * *(f4 + i);
		}
	    }
	    for (i = 0; i < Dim; i++)
		(*p)[i] = (*p)[i] * 2;
	    if (no1 == 2)
		return (FLOAT *)values;
	    p++;
	case 2:
	    a = e->verts[0];
	    b = e->verts[1];
	    c = e->verts[3];

	    if ((a < b) && (b < c)) {
		g1 = J;
		g2 = J + 4;
		g4 = J + 12;
		*f1 = *(g2 + 1) * *(g4 + 2) - *(g2 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g4) - *(g2) * *(g4 + 2);
		*(f1 + 2) = *(g2) * *(g4 + 1) - *(g2 + 1) * *(g4);
		*f2 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f2 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f4 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L0 * *(f1 + i) + L1 * *(f2 + i) + L3 * *(f4 + i);
		}
	    }
	    if ((a < c) && (c < b)) {
		g1 = J;
		g2 = J + 12;
		g4 = J + 4;
		*f1 = *(g2 + 1) * *(g4 + 2) - *(g2 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g4) - *(g2) * *(g4 + 2);
		*(f1 + 2) = *(g2) * *(g4 + 1) - *(g2 + 1) * *(g4);
		*f2 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f2 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f4 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L0 * *(f1 + i) + L3 * *(f2 + i) + L1 * *(f4 + i);
		}
	    }
	    if ((b < c) && (c < a)) {
		g1 = J + 4;
		g2 = J + 12;
		g4 = J;
		*f1 = *(g2 + 1) * *(g4 + 2) - *(g2 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g4) - *(g2) * *(g4 + 2);
		*(f1 + 2) = *(g2) * *(g4 + 1) - *(g2 + 1) * *(g4);
		*f2 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f2 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f4 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L1 * *(f1 + i) + L3 * *(f2 + i) + L0 * *(f4 + i);
		}
	    }
	    if ((b < a) && (a < c)) {
		g1 = J + 4;
		g2 = J;
		g4 = J + 12;
		*f1 = *(g2 + 1) * *(g4 + 2) - *(g2 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g4) - *(g2) * *(g4 + 2);
		*(f1 + 2) = *(g2) * *(g4 + 1) - *(g2 + 1) * *(g4);
		*f2 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f2 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f4 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L1 * *(f1 + i) + L0 * *(f2 + i) + L3 * *(f4 + i);
		}
	    }
	    if ((c < b) && (b < a)) {
		g1 = J + 12;
		g2 = J + 4;
		g4 = J;
		*f1 = *(g2 + 1) * *(g4 + 2) - *(g2 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g4) - *(g2) * *(g4 + 2);
		*(f1 + 2) = *(g2) * *(g4 + 1) - *(g2 + 1) * *(g4);
		*f2 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f2 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f4 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L3 * *(f1 + i) + L1 * *(f2 + i) + L0 * *(f4 + i);
		}
	    }
	    if ((c < a) && (a < b)) {
		g1 = J + 12;
		g2 = J;
		g4 = J + 4;
		*f1 = *(g2 + 1) * *(g4 + 2) - *(g2 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g4) - *(g2) * *(g4 + 2);
		*(f1 + 2) = *(g2) * *(g4 + 1) - *(g2 + 1) * *(g4);
		*f2 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f2 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f4 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L3 * *(f1 + i) + L0 * *(f2 + i) + L1 * *(f4 + i);
		}
	    }
	    for (i = 0; i < Dim; i++)
		(*p)[i] = (*p)[i] * 2;
	    if (no1 == 3)
		return (FLOAT *)values;
	    p++;
	case 3:
	    a = e->verts[0];
	    b = e->verts[1];
	    c = e->verts[2];

	    if ((a < b) && (b < c)) {
		g1 = J;
		g2 = J + 4;
		g3 = J + 8;
		*f1 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f1 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);
		*f2 = *(g3 + 1) * *(g1 + 2) - *(g3 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g1) - *(g3) * *(g1 + 2);
		*(f2 + 2) = *(g3) * *(g1 + 1) - *(g3 + 1) * *(g1);
		*f3 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f3 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L0 * *(f1 + i) + L1 * *(f2 + i) + L2 * *(f3 + i);
		}
	    }
	    if ((a < c) && (c < b)) {
		g1 = J;
		g2 = J + 8;
		g3 = J + 4;
		*f1 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f1 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);
		*f2 = *(g3 + 1) * *(g1 + 2) - *(g3 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g1) - *(g3) * *(g1 + 2);
		*(f2 + 2) = *(g3) * *(g1 + 1) - *(g3 + 1) * *(g1);
		*f3 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f3 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L0 * *(f1 + i) + L2 * *(f2 + i) + L1 * *(f3 + i);
		}
	    }
	    if ((b < a) && (a < c)) {
		g1 = J + 4;
		g2 = J;
		g3 = J + 8;
		*f1 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f1 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);
		*f2 = *(g3 + 1) * *(g1 + 2) - *(g3 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g1) - *(g3) * *(g1 + 2);
		*(f2 + 2) = *(g3) * *(g1 + 1) - *(g3 + 1) * *(g1);
		*f3 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f3 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L1 * *(f1 + i) + L0 * *(f2 + i) + L2 * *(f3 + i);
		}
	    }
	    if ((b < c) && (c < a)) {
		g1 = J + 4;
		g2 = J + 8;
		g3 = J;
		*f1 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f1 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);
		*f2 = *(g3 + 1) * *(g1 + 2) - *(g3 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g1) - *(g3) * *(g1 + 2);
		*(f2 + 2) = *(g3) * *(g1 + 1) - *(g3 + 1) * *(g1);
		*f3 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f3 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L1 * *(f1 + i) + L2 * *(f2 + i) + L0 * *(f3 + i);
		}
	    }
	    if ((c < a) && (a < b)) {
		g1 = J + 8;
		g2 = J;
		g3 = J + 4;
		*f1 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f1 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);
		*f2 = *(g3 + 1) * *(g1 + 2) - *(g3 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g1) - *(g3) * *(g1 + 2);
		*(f2 + 2) = *(g3) * *(g1 + 1) - *(g3 + 1) * *(g1);
		*f3 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f3 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L2 * *(f1 + i) + L0 * *(f2 + i) + L1 * *(f3 + i);
		}
	    }
	    if ((c < b) && (b < a)) {
		g1 = J + 8;
		g2 = J + 4;
		g3 = J;
		*f1 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f1 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);
		*f2 = *(g3 + 1) * *(g1 + 2) - *(g3 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g1) - *(g3) * *(g1 + 2);
		*(f2 + 2) = *(g3) * *(g1 + 1) - *(g3 + 1) * *(g1);
		*f3 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f3 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i] =
			L2 * *(f1 + i) + L1 * *(f2 + i) + L0 * *(f3 + i);
		}
	    }
	    for (i = 0; i < Dim; i++)
		(*p)[i] = (*p)[i] * 2;
	    if (no1 == 4)
		return (FLOAT *)values;
    }

//      phgFree(f1);
//      phgFree(f2);
//      phgFree(f3);
//      phgFree(f4);
    return (FLOAT *)values;
}

static const FLOAT *
RT1_grad(DOF *dof, SIMPLEX *e, int no0, int no1, const FLOAT *lambda)
{
    static FLOAT values[NFace][Dim][Dim + 1];
    FLOAT (*p)[Dim][Dim + 1] = values;
    const FLOAT *J, *g1, *g2, *g3, *g4;
//      FLOAT *f1,*f2,*f3,*f4;
    INT i, j, a, b, c;
    FLOAT f1[3], f2[3], f3[3], f4[3];
    J = phgGeomGetJacobian(dof->g, e);
//      f1 = phgAlloc(3*sizeof(FLOAT));
//      f2 = phgAlloc(3*sizeof(FLOAT));
//      f3 = phgAlloc(3*sizeof(FLOAT));
//      f4 = phgAlloc(3*sizeof(FLOAT));

    if (no1 <= 0)
	no1 = NFace;
    assert(no0 < no1);

    switch (no0) {
	case 0:

	    a = e->verts[1];
	    b = e->verts[2];
	    c = e->verts[3];
	    if ((a < b) && (b < c)) {
		g2 = J + 4;
		g3 = J + 8;
		g4 = J + 12;
		*f2 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f2 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g2 + 2) - *(g4 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g2) - *(g4) * *(g2 + 2);
		*(f3 + 2) = *(g4) * *(g2 + 1) - *(g4 + 1) * *(g2);
		*f4 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f4 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);

		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = 0;
		    (*p)[i][1] = *(f2 + i);
		    (*p)[i][2] = *(f3 + i);
		    (*p)[i][3] = *(f4 + i);
		}
	    }
	    if ((a < c) && (c < b)) {
		g2 = J + 4;
		g3 = J + 12;
		g4 = J + 8;
		*f2 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f2 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g2 + 2) - *(g4 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g2) - *(g4) * *(g2 + 2);
		*(f3 + 2) = *(g4) * *(g2 + 1) - *(g4 + 1) * *(g2);
		*f4 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f4 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);

		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = 0;
		    (*p)[i][1] = *(f2 + i);
		    (*p)[i][2] = *(f4 + i);
		    (*p)[i][3] = *(f3 + i);
		}
	    }
	    if ((b < c) && (c < a)) {
		g2 = J + 8;
		g3 = J + 12;
		g4 = J + 4;
		*f2 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f2 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g2 + 2) - *(g4 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g2) - *(g4) * *(g2 + 2);
		*(f3 + 2) = *(g4) * *(g2 + 1) - *(g4 + 1) * *(g2);
		*f4 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f4 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);

		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = 0;
		    (*p)[i][1] = *(f4 + i);
		    (*p)[i][2] = *(f2 + i);
		    (*p)[i][3] = *(f3 + i);
		}
	    }
	    if ((b < a) && (a < c)) {
		g2 = J + 8;
		g3 = J + 4;
		g4 = J + 12;
		*f2 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f2 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g2 + 2) - *(g4 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g2) - *(g4) * *(g2 + 2);
		*(f3 + 2) = *(g4) * *(g2 + 1) - *(g4 + 1) * *(g2);
		*f4 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f4 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);

		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = 0;
		    (*p)[i][1] = *(f3 + i);
		    (*p)[i][2] = *(f2 + i);
		    (*p)[i][3] = *(f4 + i);
		}
	    }
	    if ((c < a) && (a < b)) {
		g2 = J + 12;
		g3 = J + 4;
		g4 = J + 8;
		*f2 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f2 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g2 + 2) - *(g4 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g2) - *(g4) * *(g2 + 2);
		*(f3 + 2) = *(g4) * *(g2 + 1) - *(g4 + 1) * *(g2);
		*f4 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f4 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);

		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = 0;
		    (*p)[i][1] = *(f3 + i);
		    (*p)[i][2] = *(f4 + i);
		    (*p)[i][3] = *(f2 + i);
		}
	    }
	    if ((c < b) && (b < a)) {
		g2 = J + 12;
		g3 = J + 8;
		g4 = J + 4;
		*f2 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f2 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g2 + 2) - *(g4 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g2) - *(g4) * *(g2 + 2);
		*(f3 + 2) = *(g4) * *(g2 + 1) - *(g4 + 1) * *(g2);
		*f4 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f4 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);

		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = 0;
		    (*p)[i][1] = *(f4 + i);
		    (*p)[i][2] = *(f3 + i);
		    (*p)[i][3] = *(f2 + i);
		}
	    }
	    for (i = 0; i < Dim; i++)
		for (j = 0; j < Dim + 1; j++)
		    (*p)[i][j] = (*p)[i][j] * 2;
	    if (no1 == 1)
		return (FLOAT *)values;
	    p++;
	case 1:
	    a = e->verts[0];
	    b = e->verts[2];
	    c = e->verts[3];
	    if ((a < b) && (b < c)) {
		g1 = J;
		g3 = J + 8;
		g4 = J + 12;
		*f1 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f1 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f3 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g3 + 2) - *(g1 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g3) - *(g1) * *(g3 + 2);
		*(f4 + 2) = *(g1) * *(g3 + 1) - *(g1 + 1) * *(g3);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f1 + i);
		    (*p)[i][1] = 0;
		    (*p)[i][2] = *(f3 + i);
		    (*p)[i][3] = *(f4 + i);
		}
	    }
	    if ((a < c) && (c < b)) {
		g1 = J;
		g3 = J + 12;
		g4 = J + 8;
		*f1 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f1 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f3 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g3 + 2) - *(g1 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g3) - *(g1) * *(g3 + 2);
		*(f4 + 2) = *(g1) * *(g3 + 1) - *(g1 + 1) * *(g3);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f1 + i);
		    (*p)[i][1] = 0;
		    (*p)[i][2] = *(f4 + i);
		    (*p)[i][3] = *(f3 + i);
		}
	    }
	    if ((b < c) && (c < a)) {
		g1 = J + 8;
		g3 = J + 12;
		g4 = J;
		*f1 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f1 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f3 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g3 + 2) - *(g1 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g3) - *(g1) * *(g3 + 2);
		*(f4 + 2) = *(g1) * *(g3 + 1) - *(g1 + 1) * *(g3);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f4 + i);
		    (*p)[i][1] = 0;
		    (*p)[i][2] = *(f1 + i);
		    (*p)[i][3] = *(f3 + i);
		}
	    }
	    if ((b < a) && (a < c)) {
		g1 = J + 8;
		g3 = J;
		g4 = J + 12;
		*f1 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f1 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f3 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g3 + 2) - *(g1 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g3) - *(g1) * *(g3 + 2);
		*(f4 + 2) = *(g1) * *(g3 + 1) - *(g1 + 1) * *(g3);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f3 + i);
		    (*p)[i][1] = 0;
		    (*p)[i][2] = *(f1 + i);
		    (*p)[i][3] = *(f4 + i);
		}
	    }
	    if ((c < b) && (b < a)) {
		g1 = J + 12;
		g3 = J + 8;
		g4 = J;
		*f1 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f1 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f3 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g3 + 2) - *(g1 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g3) - *(g1) * *(g3 + 2);
		*(f4 + 2) = *(g1) * *(g3 + 1) - *(g1 + 1) * *(g3);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f4 + i);
		    (*p)[i][1] = 0;
		    (*p)[i][2] = *(f3 + i);
		    (*p)[i][3] = *(f1 + i);
		}
	    }
	    if ((c < a) && (a < b)) {
		g1 = J + 12;
		g3 = J;
		g4 = J + 8;
		*f1 = *(g3 + 1) * *(g4 + 2) - *(g3 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g3 + 2) * *(g4) - *(g3) * *(g4 + 2);
		*(f1 + 2) = *(g3) * *(g4 + 1) - *(g3 + 1) * *(g4);
		*f3 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f3 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f3 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g3 + 2) - *(g1 + 2) * *(g3 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g3) - *(g1) * *(g3 + 2);
		*(f4 + 2) = *(g1) * *(g3 + 1) - *(g1 + 1) * *(g3);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f3 + i);
		    (*p)[i][1] = 0;
		    (*p)[i][2] = *(f4 + i);
		    (*p)[i][3] = *(f1 + i);
		}
	    }
	    for (i = 0; i < Dim; i++)
		for (j = 0; j < Dim + 1; j++)
		    (*p)[i][j] = (*p)[i][j] * 2;
	    if (no1 == 2)
		return (FLOAT *)values;
	    p++;
	case 2:

	    a = e->verts[0];
	    b = e->verts[1];
	    c = e->verts[3];

	    if ((a < b) && (b < c)) {
		g1 = J;
		g2 = J + 4;
		g4 = J + 12;
		*f1 = *(g2 + 1) * *(g4 + 2) - *(g2 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g4) - *(g2) * *(g4 + 2);
		*(f1 + 2) = *(g2) * *(g4 + 1) - *(g2 + 1) * *(g4);
		*f2 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f2 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f4 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f1 + i);
		    (*p)[i][1] = *(f2 + i);
		    (*p)[i][2] = 0;
		    (*p)[i][3] = *(f4 + i);
		}
	    }
	    if ((a < c) && (c < b)) {
		g1 = J;
		g2 = J + 12;
		g4 = J + 4;
		*f1 = *(g2 + 1) * *(g4 + 2) - *(g2 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g4) - *(g2) * *(g4 + 2);
		*(f1 + 2) = *(g2) * *(g4 + 1) - *(g2 + 1) * *(g4);
		*f2 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f2 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f4 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f1 + i);
		    (*p)[i][1] = *(f4 + i);
		    (*p)[i][2] = 0;
		    (*p)[i][3] = *(f2 + i);
		}
	    }
	    if ((b < c) && (c < a)) {
		g1 = J + 4;
		g2 = J + 12;
		g4 = J;
		*f1 = *(g2 + 1) * *(g4 + 2) - *(g2 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g4) - *(g2) * *(g4 + 2);
		*(f1 + 2) = *(g2) * *(g4 + 1) - *(g2 + 1) * *(g4);
		*f2 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f2 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f4 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f4 + i);
		    (*p)[i][1] = *(f1 + i);
		    (*p)[i][2] = 0;
		    (*p)[i][3] = *(f2 + i);
		}
	    }
	    if ((b < a) && (a < c)) {
		g1 = J + 4;
		g2 = J;
		g4 = J + 12;
		*f1 = *(g2 + 1) * *(g4 + 2) - *(g2 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g4) - *(g2) * *(g4 + 2);
		*(f1 + 2) = *(g2) * *(g4 + 1) - *(g2 + 1) * *(g4);
		*f2 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f2 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f4 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f2 + i);
		    (*p)[i][1] = *(f1 + i);
		    (*p)[i][2] = 0;
		    (*p)[i][3] = *(f4 + i);
		}
	    }
	    if ((c < b) && (b < a)) {
		g1 = J + 12;
		g2 = J + 4;
		g4 = J;
		*f1 = *(g2 + 1) * *(g4 + 2) - *(g2 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g4) - *(g2) * *(g4 + 2);
		*(f1 + 2) = *(g2) * *(g4 + 1) - *(g2 + 1) * *(g4);
		*f2 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f2 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f4 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f4 + i);
		    (*p)[i][1] = *(f2 + i);
		    (*p)[i][2] = 0;
		    (*p)[i][3] = *(f1 + i);
		}
	    }
	    if ((c < a) && (a < b)) {
		g1 = J + 12;
		g2 = J;
		g4 = J + 4;
		*f1 = *(g2 + 1) * *(g4 + 2) - *(g2 + 2) * *(g4 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g4) - *(g2) * *(g4 + 2);
		*(f1 + 2) = *(g2) * *(g4 + 1) - *(g2 + 1) * *(g4);
		*f2 = *(g4 + 1) * *(g1 + 2) - *(g4 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g4 + 2) * *(g1) - *(g4) * *(g1 + 2);
		*(f2 + 2) = *(g4) * *(g1 + 1) - *(g4 + 1) * *(g1);
		*f4 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f4 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f4 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f2 + i);
		    (*p)[i][1] = *(f4 + i);
		    (*p)[i][2] = 0;
		    (*p)[i][3] = *(f1 + i);
		}
	    }
	    for (i = 0; i < Dim; i++)
		for (j = 0; j < Dim + 1; j++)
		    (*p)[i][j] = (*p)[i][j] * 2;
	    if (no1 == 3)
		return (FLOAT *)values;
	    p++;
	case 3:
	    a = e->verts[0];
	    b = e->verts[1];
	    c = e->verts[2];

	    if ((a < b) && (b < c)) {
		g1 = J;
		g2 = J + 4;
		g3 = J + 8;
		*f1 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f1 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);
		*f2 = *(g3 + 1) * *(g1 + 2) - *(g3 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g1) - *(g3) * *(g1 + 2);
		*(f2 + 2) = *(g3) * *(g1 + 1) - *(g3 + 1) * *(g1);
		*f3 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f3 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f1 + i);
		    (*p)[i][1] = *(f2 + i);
		    (*p)[i][2] = *(f3 + i);
		    (*p)[i][3] = 0;
		}
	    }
	    if ((a < c) && (c < b)) {
		g1 = J;
		g2 = J + 8;
		g3 = J + 4;
		*f1 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f1 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);
		*f2 = *(g3 + 1) * *(g1 + 2) - *(g3 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g1) - *(g3) * *(g1 + 2);
		*(f2 + 2) = *(g3) * *(g1 + 1) - *(g3 + 1) * *(g1);
		*f3 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f3 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);

		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f1 + i);
		    (*p)[i][1] = *(f3 + i);
		    (*p)[i][2] = *(f2 + i);
		    (*p)[i][3] = 0;
		}
	    }
	    if ((b < a) && (a < c)) {
		g1 = J + 4;
		g2 = J;
		g3 = J + 8;
		*f1 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f1 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);
		*f2 = *(g3 + 1) * *(g1 + 2) - *(g3 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g1) - *(g3) * *(g1 + 2);
		*(f2 + 2) = *(g3) * *(g1 + 1) - *(g3 + 1) * *(g1);
		*f3 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f3 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f2 + i);
		    (*p)[i][1] = *(f1 + i);
		    (*p)[i][2] = *(f3 + i);
		    (*p)[i][3] = 0;
		}
	    }
	    if ((b < c) && (c < a)) {
		g1 = J + 4;
		g2 = J + 8;
		g3 = J;
		*f1 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f1 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);
		*f2 = *(g3 + 1) * *(g1 + 2) - *(g3 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g1) - *(g3) * *(g1 + 2);
		*(f2 + 2) = *(g3) * *(g1 + 1) - *(g3 + 1) * *(g1);
		*f3 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f3 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f3 + i);
		    (*p)[i][1] = *(f1 + i);
		    (*p)[i][2] = *(f2 + i);
		    (*p)[i][3] = 0;
		}
	    }
	    if ((c < a) && (a < b)) {
		g1 = J + 8;
		g2 = J;
		g3 = J + 4;
		*f1 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f1 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);
		*f2 = *(g3 + 1) * *(g1 + 2) - *(g3 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g1) - *(g3) * *(g1 + 2);
		*(f2 + 2) = *(g3) * *(g1 + 1) - *(g3 + 1) * *(g1);
		*f3 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f3 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f2 + i);
		    (*p)[i][1] = *(f3 + i);
		    (*p)[i][2] = *(f1 + i);
		    (*p)[i][3] = 0;
		}
	    }
	    if ((c < b) && (b < a)) {
		g1 = J + 8;
		g2 = J + 4;
		g3 = J;
		*f1 = *(g2 + 1) * *(g3 + 2) - *(g2 + 2) * *(g3 + 1);
		*(f1 + 1) = *(g2 + 2) * *(g3) - *(g2) * *(g3 + 2);
		*(f1 + 2) = *(g2) * *(g3 + 1) - *(g2 + 1) * *(g3);
		*f2 = *(g3 + 1) * *(g1 + 2) - *(g3 + 2) * *(g1 + 1);
		*(f2 + 1) = *(g3 + 2) * *(g1) - *(g3) * *(g1 + 2);
		*(f2 + 2) = *(g3) * *(g1 + 1) - *(g3 + 1) * *(g1);
		*f3 = *(g1 + 1) * *(g2 + 2) - *(g1 + 2) * *(g2 + 1);
		*(f3 + 1) = *(g1 + 2) * *(g2) - *(g1) * *(g2 + 2);
		*(f3 + 2) = *(g1) * *(g2 + 1) - *(g1 + 1) * *(g2);
		for (i = 0; i < Dim; i++) {
		    (*p)[i][0] = *(f3 + i);
		    (*p)[i][1] = *(f2 + i);
		    (*p)[i][2] = *(f1 + i);
		    (*p)[i][3] = 0;
		}
	    }
	    for (i = 0; i < Dim; i++)
		for (j = 0; j < Dim + 1; j++)
		    (*p)[i][j] = (*p)[i][j] * 2;
	    if (no1 == 4)
		return (FLOAT *)values;
    }
//      phgFree(f1);
//      phgFree(f2);
//      phgFree(f3);
//      phgFree(f4);

    return (FLOAT *)values;
}

static const FLOAT *
face_normal(GRID *g, SIMPLEX *e, int face)
{
    static FLOAT normal[3 * NFace];
    FLOAT *n = normal + 3 * face;
    FLOAT x0, y0, z0, x1, y1, z1, x2, y2, z2;
    COORD *c;
    int u0, u1, u2;
    if (IsLeaf(e)) {
	memcpy(n, phgGeomGetFaceNormal(g, e, face), 3 * sizeof(*n));
    }
    else {
	GetFaceVertices(e, face, u0, u1, u2);
	x0 = (*(c = g->verts + e->verts[u0]))[0];
	y0 = (*c)[1];
	z0 = (*c)[2];
	x1 = (*(c = g->verts + e->verts[u1]))[0] - x0;
	y1 = (*c)[1] - y0;
	z1 = (*c)[2] - z0;
	x2 = (*(c = g->verts + e->verts[u2]))[0] - x0;
	y2 = (*c)[1] - y0;
	z2 = (*c)[2] - z0;
	n[0] = y1 * z2 - y2 * z1;
	n[1] = z1 * x2 - z2 * x1;
	n[2] = x1 * y2 - x2 * y1;
	/* FIXME: normalize n? */
    }
    /* FIXME: adjust direction of n? */
    return n;
}

void
RT1_init(DOF *dof, SIMPLEX *e, GTYPE type, int index,
	 DOF_USER_FUNC userfunc, DOF_USER_FUNC_LAMBDA userfunc_lambda,
	 const FLOAT *funcvalues, FLOAT *dofvalues, FLOAT **pdofvalues)
{
    GRID *g = dof->g;
    const FLOAT *n = NULL;
    const FLOAT *m = NULL;
    const FLOAT *p, *w, *bas;
    FLOAT d, d0;
    COORD *c;
    FLOAT x, y, z, *pd;
    FLOAT lambda[Dim + 1];
    INT i, j, k = 0, v0, v1, v2, face, order, nn;
    DOF_TYPE *t;
    QUAD *quad;
    FLOAT values[Dim], A[NFace][NFace], b[NFace];

    i = DofTypeOrder(dof, e);
    order = i + i;
    quad = phgQuadGetQuad2D(order);
    Unused(pdofvalues);

    switch (type) {
	case FACE:
	    t = dof->type;
	    nn = t->np_face;
	    break;
	default:
	    t = NULL;
	    nn = 0;
    }

    if (funcvalues != NULL) {
	memcpy(dofvalues, funcvalues, nn * dof->dim * sizeof(*dofvalues));
	return;
    }

    if (userfunc == NULL && userfunc_lambda == NULL)
	phgError(1, "%s:%d, invalid user function!\n", __FILE__, __LINE__);

    switch (type) {
	case FACE:
	    i = DofTypeOrder(dof, e);
	    order = i + i;
	    quad = phgQuadGetQuad2D(order);

	    n = face_normal(g, e, index);
	    v0 = GetFaceVertex(index, 0);
	    v1 = GetFaceVertex(index, 1);
	    v2 = GetFaceVertex(index, 2);
	    lambda[index] = 0.;
	    p = quad->points;
	    w = quad->weights;
	    d = 0.;
	    for (i = 0; i < quad->npoints; i++) {
		m = n;
		lambda[v0] = *(p++);
		lambda[v1] = *(p++);
		lambda[v2] = *(p++);
		if (userfunc_lambda != NULL) {
		    userfunc_lambda(dof, e, index, lambda, values);
		}
		else {
		    x = y = z = 0.;
		    for (k = 0; k < NVert; k++) {
			c = g->verts + e->verts[k];
			x += (*c)[0] * lambda[k];
			y += (*c)[1] * lambda[k];
			z += (*c)[2] * lambda[k];
		    }
		    userfunc(x, y, z, values);
		}
		pd = values;
		d0 = 0.;
		for (k = 0; k < Dim; k++) {
		    d0 += *(pd++) * *(m++);
		}
		d += d0 * *(w++);
	    }
	    b[index] = d;
	    break;
	default:
	    break;
    }
    switch (type) {
	case FACE:
	    *dofvalues = b[index] * phgGeomGetFaceArea(g, e, index);
	    break;
	default:
	    break;
    }

}

static BYTE RT1_orders[] = { 1, 1, 1, 1 };	/* Right??? */

static DOF_TYPE DOF_RT1_ =
    { DofCache, "RT1", NULL, RT1_orders, DOF_DG0, NULL, NULL,
    phgDofInterpC2FGeneric, phgDofInterpF2CGeneric, RT1_init, RT1_bas,
	RT1_grad,
    NULL, FE_Hdiv,
    FALSE, FALSE, -1,		/* invariant, free_after_use, id */
    NFace, 1, -1, 3,		/* nbas, order, cont, dim */
    0, 0, 1, 0
};				/* np_vert, np_edge, np_face, np_elem */

